#ifndef READ_DATA_H
#define READ_DATA_H
//#include "cg_coeff.h"
//#include "const_nucl.h"
//#include "variables.h" 
#include "../DF_Common/basis_state.h"
#include "../DF_Common/com_sparse_row.h"
#include "../DF_Common/coupling_coeff.h"
#include <string>
#include <vector>

using std::cout;
using std::endl;
using std::string;

class SinglP_RD : public SingleP {
public:
  double hbar_omega;
  int nshell;
  int orbits_num;
  vector<int> rd2out;
  vector<int> out2rd;
  void build_from_Motern(string filename);
  void match_sp(vector<State_SP> state_out);
  void sort_Luigi();
};

struct TBME_Conf {
  int F;
  int I;
  double val;
  int flag;
};

class DoubleP_RD : public DoubleP {
public:
  vector<vector<TBME_Conf>> conf; //[ch]{F,I,val,flag}
  SMatrix_Ch conf_vec_SM;         //[ch][F][I] = conf_index
  int Conf_Num;
  void build_for_Motern(SinglP_RD &sp_t);
  void build_RD_configure_Mortern();
  void build_for_Luigi_T(SinglP_RD &sp_t);
  void build_RD_configure_Luigi();
};

class Read_Data_II {
public:
  SinglP_RD sp_rd;
  DoubleP_RD dp_rd;
  CG_coeff cg_rd;
  int Conf_Num;
  double EPS = 0.00000000001;
  void read_orbits(string file);
  void read_int(string file);

  void read_orbits_Luigi(int E_sp);
  void read_int_Luigi(string file);
  void check_int_Luigi(string file_pn, string file_out);

  double get_tbme_T(int &a, int &b, int &c, int &d, int &J2, int &T2);
  double get_tbme_Tz_from_T(int &a, int &b, int &c, int &d, int &J2);
  double get_tbme_Tz(int &a, int &b, int &c, int &d, int &J2);
  double get_tbme(int &a, int &b, int &c, int &d, int &J2) {
    // double res = get_tbme_Tz_from_T(a, b, c, d, J2);
    double res = get_tbme_Tz(a, b, c, d, J2);
    return res;
  }
  double get_tbme_unnormalized(int &a, int &b, int &c, int &d, int &J2) {
    double res = get_tbme_unnormalized_T(a, b, c, d, J2);
    return res;
  }
  double get_tbme_unnormalized_Tz(int &a, int &b, int &c, int &d, int &J2) {
    double tbme = this->get_tbme(a, b, c, d, J2);
    double Norm_fac = 1.0;
    double Norm_fac_t = 1.0;
    if (a == b)
      Norm_fac_t = Norm_fac_t * (std::sqrt(2));
    if (c == d)
      Norm_fac_t = Norm_fac_t * (std::sqrt(2));
    Norm_fac = 1.0 / (Norm_fac_t);
    double tbme_un = tbme / Norm_fac;
    return tbme_un;
  }
  double get_tbme_unnormalized_T(int &a, int &b, int &c, int &d, int &J2) {
    double tbme = this->get_tbme(a, b, c, d, J2);
    double Norm_fac = 1.0;
    double Norm_fac_t = 1.0;
    if (same_group(a - 1, b - 1) == 1)
      Norm_fac_t = Norm_fac_t * (std::sqrt(2));
    if (same_group(c - 1, d - 1) == 1)
      Norm_fac_t = Norm_fac_t * (std::sqrt(2));
    Norm_fac = 1.0 / (Norm_fac_t);
    double tbme_un = tbme / Norm_fac;
    return tbme_un;
  }
  double get_tbme_this_T(int &a, int &b, int &c, int &d, int &ch);

private:
  double get_tbme_this(int &a, int &b, int &c, int &d, int &ch);

  int phase(int n) { return n % 2 ? -1 : 1; }
  int same_group(int a_tz, int b_tz);
};

struct state {
  int num_j;
  int num_m;
  int nr;
  int l;
  int j2;
  int tz;
  int m2;
  int P_H; // P_H= 0 particle, 1 hole state
  double energy;
};

struct twoBodyelement {
  int num;
  int Tz;
  int Par;
  int J2;
  int a, b, c, d;
  double V, Hcom, rirj, pipj;
};
/**/
struct twoBody_index_ab {
  int J2, M2;
  int Tz;
  int a, b;
  int index;
};

class Read_Data {
public:
  // Read_Data(std::string orbits_file,std::string interaction_file)
  void ReadOrbits(std::string file_name);
  void ReadInteraction(std::string file_name);
  double get2Body_j(int &a, int &b, int &c, int &d, int &J2,
                    int SphericalHf_info, int hw_A);
  double get2Body_j_unnormalized(int &a, int &b, int &c, int &d, int &J2,
                                 int SphericalHf_info, int hw_A) {
    double tbme = this->get2Body_j(a, b, c, d, J2, SphericalHf_info, hw_A);
    double Norm_fac = 1.0;
    if (a == b)
      Norm_fac = Norm_fac / (std::sqrt(2));
    if (c == d)
      Norm_fac = Norm_fac / (std::sqrt(2));
    double tbme_un = tbme / Norm_fac;
    return tbme_un;
  }
  void check(std::string file_name, std::string file_name2);
  // double get2Body(int &a,int &b,int &c,int &d,int &J2,double &g2Bv,int
  // &Tz,int &Par);
  int nshell_r, h_w_r, orbits_num;
  int J2max, num_2Be;

  CG_coeff cg_quick; // initial in ReadOrbits
  // cg_quick.init_cal(J2max);
  std::vector<state> J_basis;
  std::vector<twoBodyelement> data_v;

  // for m-scheme !!
  void init_M_basis();
  double get2Body_m(int &a, int &b, int &c, int &d);
  std::vector<state> M_basis;
  std::vector<std::vector<int>> tBiab_M_mat;
  //
  std::vector<std::vector<std::vector<std::vector<double>>>> m_tBiV;

private:
  // this variables are used to find 2-body-element quickly
  // for example:
  //  index1=tBiab_mat[j2_2][a-1][b-1];
  //  index2=tBiab_mat[j2_2][c-1][d-1];
  //  v_index=tBiV[j2_2][index1][index2];
  //  g2Bv=tBe_J[j2_2][v_index].V-tBe_J[j2_2][v_index].pipj*hw_A;

  // classify data_v by J2
  std::vector<std::vector<twoBodyelement>> tBe_J;
  //
  std::vector<std::vector<std::vector<int>>> tBiab_mat;
  //
  std::vector<std::vector<std::vector<int>>> tBiV;

  double get2Body(int &a, int &b, int &c, int &d, int &J2, double &g2Bv,
                  int &Tz, int &Par, int hw_A);

  // for   M scheme !!!
  // this variables are used to find 2-body-element quickly in M scheme
  // std::vector<std::vector<int> > tBiab_M_mat;
  //
  // std::vector<std::vector<std::vector<std::vector<double> > > > m_tBiV;

  double get2Body_mt(int &a, int &b, int &c, int &d, int hw_A);

  void find_m_tBiab(int &tz, int &M2,
                    std::vector<std::vector<int>> &ab_same_M2tz);
};

#endif
